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Abstract 

The accuracy of the energy landscape of silicon systems obtained from various density functional 
methods, a tight binding scheme and force fields is studied. Quantum Monte Carlo results serve as 
quasi exact reference values. In addition to the well known accuracy of DFT methods for geometric 
ground states and metastable configurations we find that DFT methods give a similar accuracy for 
transition states and thus a good overall description of the energy landscape. On the other hand, 
force fields give a very poor description of the landscape that are in most cases too rugged and 
contain many fake local minima and saddle points or ones that have the wrong height. 
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I. INTRODUCTION 



In spite of the great progress in density functional methods for treating large systems, 
it is at present not possible to treat systems with more than about 1000 atoms in complex 
simulations where forces and energies have to be evaluated many times. This is for instance 
necessary in molecular dynamics simulations where one has to follow the evolution of the 
system over long time intervals or in global optimization methods for finding the ground 
state geometry. In these kinds of situations faster and more approximate methods such as 
force fields or tight binding schemes are widely used. Because of its technological impor- 
tance several widely used force fields exist for silicon and large scale simulations, which are 
not feasible with density functional methods, are frequently performed using these more 
approximate methods. 

These force fields are typically fitted to a data set of ground state structures, usually 
containing crystalline structures and sometimes also non-periodic structures. An accurate 
description of some ground state geometries is however not sufficient to ensure accurate 
dynamical simulations. Dynamical properties such as diffusion coefficients are related to 
other properties of the energy landscape such as barrier heights. The distribution of bar- 
rier heights and other properties of the silicon potential energy surface (PES) have been 
studied using forcefieldsi. In this paper we study overall properties of the energy landscape 
which are relevant in many different contexts. We look in particular at the accuracy of the 
barrier heights in the various schemes used for large scale simulations of silicon systems. 
Since it is known that the barrier heights relevant to chemical reactions are not very well 
described with standard density functionals such as local density approximation (LDA)2 or 
Perdew-Burke-Ernzerhof (PBE),- we benchmark some of the energy barriers using accurate 
quantum Monte Carlo methods. There are two forms of QMC methods that are used for 
electronic structure calculations, the simpler variational Monte Carlo (VMC) and the more 
sophisticated diffusion Monte Carlo (DMC). In VMC, quantum mechanical expectation val- 
ues are calculated using Monte Carlo techniques to evaluate the many- dimensional integrals. 
The accuracy of the results depends crucially on the quality of the trial wave function. The 
DMC removes most of the error in the trial wavefunction. DMC is a stochastic projector 
method that projects out the ground state from the trial wavefunction using an integral 
form of the imaginary-time Schrodinger equation. For Fermionic systems, the antisymme- 
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try constraint leads to the Fermion-sign problem that is cured by fixing the nodes of the 
projected state to be those of an approximate trial wavefunction. The resulting fixed-node 
error is the main uncontrolled error in DMC. Currently, a systematic improvement of the 
wavefunction by optimization of an increasing number of variational parameters is the most 
practical approach for reducing the fixed-node error.-^ 

Force fields and other approximate methods are sometimes applied to systems that are 
very different from the systems that were in the fitting database. The question is therefore 
how reliable are force field based structure predictions of complex structures such as defects, 
interfaces or clusters. In most studies of such systems only one force field was used but 
in some exceptionally careful studies, such as in the study of dislocation kinks in silicon^ 
and a comparative study of silicon empirical interatomic potentials^, the results of several 
force fields were compared and significant discrepancies in the results obtained from different 
force fields were indeed found. To shed light on the accuracy of force fields we examine the 
configurational density of states obtained by various approximate schemes for silicon. Ideally 
there would be a one-to-one mapping between stable structures (local minima) obtained 
using approximate and accurate methods. Consequently the density of configurations per 
energy interval would be identical for approximate and accurate methods. 



II. METHODS 



In our study we have included the most common force fields for silicon, namely the Ter- 
soff force field,- the Stillinger- Weber force field,— the environmental- dependent interaction 
potential (EDIP) force fields and the Lenosky force field. The Tersoff force field, which has 
infinite range, was smoothly extrapolated to zero by a third order polynomial using cutoff 
radii that are large enough to ensure a smooth behavior of the potential. The cutoff values 
were 2.7 and 3.3 A where the smaller value denotes the radius where the polynomial takes 
over and the larger value the radius where the polynomial interaction drops to zero. As a 
representative for tight binding schemes we use the Lenosky tight (LTB) binding scheme.— 
The QMC calculations are performed using the CHAMP code developed by Umrigar, Fil- 
ippi and Toulouse. The Is, 2s, and 2p electrons of Si are eliminated using a Hartree-Fock 
pseudopotential.— The trial wave function consists of a sum of Slater determinants of single- 
particle orbitals multiplied by a Jastrow factor. The orbitals of the Slater determinant are 
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taken from a DFT calculation using GAMESS^ 1 ^ with the B3LYP exchange-correlation 
functional. The excitations included in the sum of Slater determinants are those with the 
largest weights in a configuration interaction with single and double excitation (CISD) calcu- 
lation. Configuration state functions (CSFs), i.e., linear combinations of determinants that 
have the correct spatial and spin symmetries, are used to reduce the number of variational 
parameters. The Jastrow correlation function describes electron-electron, electron-nuclear 
and electron-electron nuclear correlations. The Jastrow parameters and the CSF coefficients 
are optimized in variational Monte Carlo (VMC) using a recently developed energy mini- 
mization method.-^ Finally, diffusion Monte Carlo (DMC) calculations using the optimized 
trial wavefunction and a time step of 0.01 Ha -1 determine the energies of the structures. 
Most of the calculations employed a single determinant, but to estimate the size of the 
fixed-node error we performed, for some structures, VMC and DMC calculations with trial 
wavefunctions containing an increasing number of Slater determinants up to 150. Apart 
from the preparation of trial wave function for QMC calculations and DFT barrier heights 
within B3LYP, all other DFT calculations are performed with the BigDFT— package, a 
pseudopotential based^ 1 ^ DFT code with a wavelet basis set. Wavelets are a systematic 
basis set and the basis size was chosen sufficiently large that energies were converged to 
better than to 1 mHa. 

III. TRANSITION STATES 

According to transition state theory, the height and the shape of saddle points determine 
the dynamical behavior of a system.— Not much effort has been made to assess the effect of 
approximate energy landscapes on the dynamics of silicon systems. In order to investigate 
the quality of the silicon potential energy landscape within various schemes we performed 
simulations to find saddle points of small Sis clusters using DFT within the local density 
approximation. We then compared the barrier heights for these LDA configurations with 
the heights obtained from other DFT schemes, namely with a generalized gradient approxi- 
mation functional (PBE), a hybrid functional (B3LYP), and accurate quantum Monte Carlo 
(QMC) methods. The 8 atom silicon cluster was chosen because for this size, QMC is com- 
putationally not too expensive. In agreement with previous work— we found that the saddle 
point geometries are nearly identical within different DFT functionals. This justifies the use 
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(a) minimuml (b) saddle (c) minimum2 

FIG. 1: (color online) The saddle point SP2 (b) and the two neighboring minima (a), (c). It is 
obvious that the silicon atoms in the saddle point and the minima configurations are all in a similar 
environment, in particular the saddle point is very similar to the structure (c). 

of LDA geometries in all energy calculations. We also compared DFT saddle point results 
with the Lenosky tight-binding (TB) method. Quantum Monte Carlo has demonstrated its 
ability to provide accurate reaction barrier heights. Calculated values agreed with the exper- 
imental values to within the statistical error bar of 0.07 eV,— for some organic molecules and 
to within 0.005 eV (see Refs. I22|j23l ) for the well-known exchange reaction H+H2 — > H2+H. 
We will therefore consider in the following our QMC results as accurate reference values. 

To generate our saddle point configurations, we started with the putative global minimum 
of the Si§ cluster— and using the improved dimer method^ we found six saddle points (SP1 
to SP6) which connect the global minimum state to other ones or to itself - the latter 
corresponds to the exchange of two silicon atoms (SP1 and SP4). Furthermore, two more 
saddle points (SP7, SP8) are obtained starting from the first low-lying isomer. One of them 
(SP7) corresponds to the exchange of two atoms. Finding the saddle points and the adjacent 
minima was done within LDA using the BigDFT~ package. 



There are numerous publications (for a survey see Ref. l26l ) in which it is shown that 
DFT schemes give poor transition state barrier heights for chemical reactions.— Within 
LDA and GGA's the results are most unsatisfactory for hydrogen transfer reactions where 
covalent hydrogen bonds are broken and formed. The most simple and prominent example 
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(a) SP1 



(b) SP2 




(g) SP7 (h) SP8 

FIG. 2: (color online) The eight Sis saddle points obtained by LDA calculations. 
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TABLE I: Comparison of the barrier height (BH) energies of eight saddle point (SP) configurations 
relative to the two neighboring minima calculated with various DFT schemes and QMC. The saddle 
point and minima configurations are from LDA. The abbreviation for the various DFT functionals 
are defined in the text. VMC and DMC stand for variational and diffusion QMC, MD-VMC 
and MD-DMC for multi-determinant VMC and DMC. The HOMO-LUMO gaps for the minima 
(HLGM) and for the saddle points (HLGS) are calculated within local density approximations 
using the BigDFT code. The statistical errors are given in parentheses. All energies are in electron 
volts. 

system HLGM HLGS LDA PBE B3LYP VMC DMC MD-VMC MD-DMC TB 



SP1-BH1/2 1.45 0.55 0.359 0.367 0.405 0.434(23) 0.406(13) 0.329(9) 0.338(8) -0.046 

SP2-BH1 1.45 0.41 0.735 0.729 0.770 0.809(20) 0.804(18) -0.286 

SP2-BH2 0.22 0.077 0.094 0.058 0.103(22) 0.069(18) 0.166 

SP3-BH1 0.96 0.64 0.043 0.050 0.056 0.046(18) 0.028(15) 0.042(9) 0.028(8) 0.201 

SP3-BH2 1.45 2.900 2.689 2.328 2.927(17) 2.845(15) 2.969(8) 2.826(8) 1.934 

SP4-BH1/2 1.45 0.13 1.065 1.053 1.075 1.233(18) 1.181(15) 1.131(10) 1.134(7) 0.496 

SP5-BH1 1.21 0.42 0.338 0.346 0.324 0.410(18) 0.378(14) 0.398 

SP5-BH2 1.45 0.766 0.761 0.800 0.883(17) 0.862(16) 0.019 

SP6-BH1 1.49 0.69 0.581 0.573 0.514 0.503(18) 0.536(16) -0.082 

SP6-BH2 1.45 1.198 1.158 1.015 1.283(18) 1.194(16) 0.921 

SP7-BH1/2 1.12 0.58 0.289 0.272 0.192 0.228(18) 0.224(16) 0.157 

SP8-BH1 0.22 0.82 0.212 0.198 0.041 0.143(17) 0.120(16) 0.144 

SP8-BH2 1.12 0.445 0.420 0.406 0.491(18) 0.445(16) -0.264 



is the exchange reaction H + Hi — > Hi + H where the DFT schemes do not predict a 
barrier at all. The poor performance seems to be due to poor cancellation of the electrostatic 
self-interaction errors in the DFT schemes.— In the literature this problem is known as "self- 
interaction error" which is related to the derealization error.— Hybrid functionals, which 
give a better error cancellation, give improved barrier heights in this case.— Nevertheless 
researchers usually resort to wavefunction methods if highly accurate barrier heights are 
needed for chemical reactions. 

Table [I] shows that in our case the situation is entirely different. Already at the LDA level 
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the barrier heights are reasonably accurate and actually slightly better than the B3LYP and 
PBE barriers. How can this surprising accuracy be explained? In contrast to chemical reac- 
tions our clusters are never torn apart into fragments when they move along the minimum 
energy pathway from one local minimum over a saddle point into another local minimum. 
Even at the transition state (see Figs. [T] and [2]) the silicon atoms are all in an environ- 
ment that is similar to the environment at a local minimum and one cannot distinguish a 
saddle point configuration from a local minimum energy configuration by inspection. DFT 
self-interaction errors^ are therefore expected to cancel to a large degree. 

Since DFT is essentially a one determinant method one would expect that DFT results are 
particularly poor when transition states have multi-determinant character. This is frequently 
the case in chemical reactions and under such circumstances multi-reference wavefunction 
methods have to be employed if accurate barrier heights are needed. Small HOMO-LUMO 
gaps are an indication of the importance of multi-reference configurations. The HOMO- 
LUMO gaps of all saddle points and the adjacent minima are presented in Table HI The 
average HOMO-LUMO gap of saddle points is less than that of the minima by ~ 0.65 eV. As 
usual the LDA and GGA gaps are much smaller than the B3LYP gaps which are certainly 
more accurate. For the configurations with small HOMO-LUMO gaps we also did multi- 
determinant QMC calculations with as many as 150 determinants. The DMC energies went 
down by no more than 2 mHa. In Fig. HI the correlation between the HOMO-LUMO gap 
and the energy change from single- to multi-determinant calculations is illustrated for seven 
configurations, in particular for the three configurations that have a HOMO-LUMO gap of 
less than 0.25 eV. The results show that the influence of a multi determinant wavefunction on 
the barrier height is very small. Another indication that the multi-determinant character of 
the saddle point configurations can be neglected is the fact the natural occupation numbers 
drop rapidly from one to zero. The occupation numbers obtained from CISD calculations 
using GAMESS with about 70 determinants are shown in Table HT1 

Table [I] also shows that the tight binding barrier heights are not reliable. The situation 
is yet worse for the force fields and we have not even attempted to give error bars. An 
additional complication, which will be discussed in next section, is that the potential energy 
surfaces of force fields contain many fake minima and consequently also many fake saddle 
points connecting the fake minima. 
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TABLE II: The occupation number of the four highest occupied orbitals and the four lowest 



unoccupied orbitals of Sis clusters obtained from CISD calculations. 
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FIG. 3: (color online) The correlation of barrier heights from various DFT functionals with the 
DMC barrier heights. 
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FIG. 4: (color online) The correlation between the energy change from single- to multi-determinant 
diffusion Monte Carlo calculations with the HOMO-LUMO gap. The (red) circles are for saddle 
point structures and the (blue) squares are for minimum point structures. 

IV. LOW ENERGY CONFIGURATIONS OF SILICON CLUSTERS AND CRYS- 
TALLINE STRUCTURES 



A. Funnel-like structure of the PES 

The potential energy surface (PES) of the Sii6 cluster was explored systematically with all 
of the aforementioned classical many-body potentials and the Tight-Binding scheme using 
the minima hopping methocP (MHM). The MHM consists of a sequence of consecutive short 
molecular dynamics runs and geometry relaxations. The MHM is a method that determines 
the global minimum of a given PES as well as other low-lying energy configurations very 
efficiently. 

Atomic clusters range from structure seekers with well denned ground state configurations 
which can be found very rapidly with the MHM to glass-like systems for which it is very 
difficult to find the ground state. The speed with which a system finds its ground state is 
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evidently a physical property of the system and should carry over to most computational 
geometry-opt imiz at ion algor it hms . 

However, we have found considerable differences in the speed of finding the ground state 
configuration with the MHM when using the various potentials to describe the Sii6 clus- 
ter. Table UTTl gives the average number umin of minima visited before finding the global 
minimum. The differences in umin can be ascribed to the configurational density of states 
(C-DOS) of the local minima for the particular potentials, discussed in section HV El A large 
n MiN indicates a high C-DOS in the low energy region and vice versa. The Sii6 cluster for 
instance looks like a structure seeker with the Lenosky force field and more like a glassy 
system with the Tersoff force field. In agreement with the observation in Fig. M the M EAM 
(Modified Embedded Atom) ansatz^ of the Lenosky force field gives a relatively smooth 
potential energy surface. 

TABLE III: Average values in 100 MHM runs, umin is the number of minima visited before 
finding the ground state configuration. 



Method nMiN 

EDIP 85 

Lenosky 10 

Reparametrized Lenosky 8 

Tersoff 116 

Stillinger- Weber 31 

Lenosky Tight-Binding 42 

DFT 32 



B. Ground state and low energy configurations of Sii6 isomers 

A database of stable configurations was generated by visiting 1000 different local minima 
with the MHM in each potential. To verify the accuracy of the potentials the ten energeti- 
cally lowest structures were relaxed to the nearest local DFT minimum. The results of the 
investigation will be discussed separately for each potential and are represented in Fig. [HJ 
The geometrical features to characterize surface properties of Siig isomers are described in 
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(a) (b) 

FIG. 5: (color online) (a) A surface atom forms a sharp corner if it defines an acute-angled tetra- 
hedral (green) or pyramidal (red) structure together with its three or four nearest neighbors re- 
spectively, (b) Two neighboring surface atoms form a sharp edge if they make an acute-angled 
tetrahedral structure with the two common nearest neighbors (blue). Four neighboring surface 
atoms are part of a facet if they form a plane (orange). 

Fig. El 

EDIP: The ground state configuration is a hollow oblate spheroid consisting of four 
parallel planes. Each plane contains four atoms forming a square which are rotated by 
45 ° with respect to the neighboring planes. The top and bottom squares have an edge 
length 2.49A, the intermediate planes 4.27A. In this configuration all atoms are five fold 
coordinated with an average bond length of 2.50A. The first nine excited configurations 
have the same general features as the ground state and only differ by forming or breaking of 
up to three bonds. They cover an energy range of 0.5 eV. When relaxed in DFT most of the 
structures are heavily deformed. In general the void regions within the structure collapse 
and the shape has a tendency to get elongated. 7 structures show sharp corners and edges 
and two structures exhibit extended facets. Only one minimum was found to be stable in 
DFT. 

Lenosky: The third and the fourth lowest energy configurations are oblate spheroid 
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and are the same as the second lowest and the global minimum configurations with EDIP 
potential respectively. However, the ground state configuration with the Lenosky potential is 
highly spherical consisting of four hexagonal curved panels with six-fold coordinated center 
atoms. Although maintaining the general form, the relaxation in DFT reveals that three 
of the four panels are transformed into planes. Remarkably, an excited configuration was 
found to relax in DFT to the above mentioned ground state configuration (energy difference 
in DFT 0.7 eV). Furthermore, there are four spherical and two elongated hollow structures, 
only one of which is stable with DFT. All other excited states are deformed heavily and the 
majority show both sharp corners and edges, covering an energy range of 0.4 eV. 

Stillinger- Weber: In contrast to the previous two potentials most configurations are 
highly elongated and often contain pentagonal elements. There are only three exceptions in- 
cluding the ground state configuration which has a hollow elliptical shape formed by 2 square 
and 8 pentagonal planes. Furthermore, none of the geometries contains over-coordinated 
atoms. The structures of all minima found with the Stillinger- Weber potential have a strong 
tendency to contain a large number of sharp corners and only few facets when relaxed in 
DFT. Although none of the structures are stable in DFT the general elongated form is often 
conserved and leads to low DFT energies, thus indicating an accurate description of the low 
energy regions on the PES. The configurational energies are scattered over a range of 0.80 
eV. 

Tersoff: The ground state geometry of the Tersoff potential is identical to the one found 
with the Stillinger- Weber potential. Only one of the nine lowest energy configurations other 
than ground state is elongated, the other eight are very similar to the ground state with 
hollow spherical shapes. The ninth excited state is 1.4 eV above the ground state. Similar 
to the Stillinger- Weber potential the structures do not include over-coordinated atoms. All 
structures are deformed after geometry optimizations with DFT and have a large number 
of sharp corners. 

Lenosky Tight-Binding: The Lenosky Tight-Binding scheme predicts a global mini- 
mum configuration which is a slightly distorted Stillinger- Weber and Tersoff ground state 
geometry. In contrast to the classical potentials both hollow elliptical and elongated struc- 
tures without void regions are predicted in equal amounts. The energy of the ninth excited 
state lies only 0.15 eV above the ground state geometry, indicating an overall shallow PES. 
Although some bond lengths are overestimated, all structures with only one exception were 
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found to be stable in DFT calculations. However, three configurations with similar geome- 
tries converged to the same minimum structure. The ordering of the minima with respect 
to the energies within the Tight-Binding scheme and the DFT calculations is in fairly good 
agreement with the ideal correlation (see Fig. [H]). While all four classical potentials fail to 
predict stable low-lying Sii6 isomers the Lenosky Tight-Binding scheme succeeds in most 
cases. 

C. Flat regions of the PES 

During DFT geometry relaxations one can encounter cases where the cluster is distorted 
considerably even though the energy decreases only slightly. Within these flat regions the 
norm of the force is small but may increase while the monotonous downhill progress in energy 
is preserved. Many steps are necessary in the steepest descent DFT geometry relaxation to 
overcome these flat plateau regions. Only the Lenosky Tight-Binding scheme provides an 
accurate energy trend when following the DFT relaxation pathway. All classical potentials 
fail to even describe the lowering of the configurational energy along the pathway (see Fig. [9]). 




123456789 
Rank in Tight-Binding 

FIG. 6: (color online) Ordering of the local minima energies in Lenosky tight-binding and DFT. 
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(a)EDIP 



(b)Lenosky 




(c)Stillinger-Weber/Tersoff (d)Lenosky Tight-Binding 

FIG. 7: (color online) Global minimum configuration of all five potentials found with the minima 
hopping method. 

With the exception of the Lenosky force field these potentials give a strongly oscillating 
energy surface instead of a flat one along the pathway. This is a first indication that the 
classical potentials give a too rough PES. The MEAM ansatz of the Lenosky force field 
seems to give smoother surfaces than the other classical potentials. 
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TABLE IV: Statistical data related to the Hessian matrix around local minima of up to 120 random 
configurations (samples) of a S130 cluster. The second line contains the number of steps used in the 
geometry relaxation, a combination of preconditioned steepest descent and preconditioned DIIS 
method (except for the DFT). The corresponding average, minimum and maximum of the largest 
and smallest eigenvalues of the Hessian matrices are listed in the lines 3 to 8 (in eV/A 2 ). The last 
line contains the average condition number k of the corresponding Hessian matrices. *Values only 
accurate in the first decimal place. 
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Furthermore, 120 random configurations were relaxed using the different potentials and 
the largest eigenvalues of the Hessian matrix were calculated when the local minimum con- 
figuration was reached. The larger and highly scattered eigenvalues (Table HVl) found within 
the classical potentials indicate that these plateau-like energy landscapes are poorly de- 
scribed by classical force fields. The Stillinger- Weber potential provides the results closest 
to the Lenosky Tight-Binding scheme, the most accurate among the potentials. This is in 
agreement with the accurate overall description of low-lying geometries with the Stillinger- 
Weber potential. The Tersoff potential requires only a small number of relaxation steps, 
indicating a high number of local minima and the absence of plateau regions on the PES. 
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TABLE V: The results often configurations of each potential relaxed with DFT. The second column 
shows the number of stable structures. The following columns show the number of structures which 
relax to the bulk crystal, to a single FFCD or two FFCD which are either neighboring (n-FFCD) 



or distant (d-FFCD). 

Method Stable Bulk FFCD n-FFCD d-FFCD 

EDIP 2 6 4 

Lenosky 10 1 2 7 

SW 7 14 5 

Tersoff 2 4 6 

LTB 10 1 1 4 4 



D. Defects in crystalline Silicon 

The MHM was used to explore the low energy region on the PES of bulk silicon. Starting 
with crystalline cubic diamond structure consisting of 216 Si atoms, 200000 local minima 
were found successively for each classical potential during the simulation. For the Lenosky 
Tight-Binding scheme only 25000 structures could be found due to limited computer time. 
Periodic boundary conditions with respect to the ground state geometry were used to provide 
the appropriate bulk conditions. The ten energetically lowest configurations of each potential 
were used as input configurations for geometry relaxations in DFT. 

The correct ground state geometry, the well-known diamond structure, is predicted with 
all the potentials. However, the structures of the first excited state of different force fields 
do not coincide. For all potentials except the Lenosky force field it is a single four fold 
coordinate defect^ (FFCD). The Lenosky potential on the other hand predicts a pair of 
two four fold coordinated defects^ in different regions of the cell as the lowest energy defect 
structure. The double FFCD is 3.99 eV higher in energy compared to the diamond structure. 

The majority of the eight other low energy geometries in the EDIP potential are structures 
containing single displaced atoms which are either four or five fold coordinated. Similar 
structures can be found with the Tersoff potential. All of these excited configurations are 
unstable in DFT calculations. The Tersoff potential additionally has minima at a variety of 
slightly distorted FFCDs which are unstable in DFT. 
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In contrast to the other three force fields, the Lenosky force field always predicts pairs 
of FFCDs as low-lying energy configurations. They are either neighboring and share a 
common atom or are distant, i.e., located in different regions of the cell. Even though the 
single FFCD, which must be the first excited state, is not predicted by the Lenosky force 
field^, all other low-lying energy configurations from the second to ninth excited states 
are stable in DFT geometry optimization. Nonetheless, the sequence with respect to the 
energy does not coincide with the sequence in DFT energies. The St illinger- Weber force 
field behaves very similarly. Only three structures were found to be unstable, the other five 
excited states all contain two interacting FFCDs. 

The best accuracy can be found with the Lenosky Tight-Binding scheme. All structures 
exist on the DFT Born-Oppenheimer surface and the energy sequence is correctly described 
with the exception of the 9 th and 10 th excited states. They are exchanged in sequence and 
show an energy difference of 0.02 eV with the Lenosky Tight-Binding scheme and 0.04 eV 
when calculated with DFT. 

E. Configurational density of states of local minima 

To describe the overall characteristics of the potential energy surface we chose the config- 
urational density of states (C-DOS), i.e., the number of configurations per energy interval. 
We approximate this C-DOS by the minima hopping density of states (MH-DOS) which is 
obtained simply by sampling the low energy region with the MHM and counting the number 
of distinct minima found in an energy interval. It has to be stressed that, in the plots we 
present in this paper, a more or less complete sampling of all minima can only be achieved 
in a very small interval around the global minimum. Only in this small interval of several 
eV we observe in the MH-DOS the expected exponential growth of the number of local min- 
ima with respect to the energy of the C-DOS. In our plots we show however a much larger 
energy interval where the number of states is the true number of states multiplied by the 
probability that a configuration in this energy range will be visited. Since this probability 
decreases with increasing energy the MH-DOS tends to zero for large energies in all our 
plots whereas the C-DOS would be orders of magnitude larger. Since the minima hopping 
method maps out higher and higher energy configurations when the minima hopping run 
is allowed to continue longer and longer, we can can compare the results of our standard 
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length MHM run, where 200 000 configurations were found, to the results of a longer run, 
where 1000000 configurations were found. As one can see from Fig. [10] the MH-DOS and 
C-DOS agree only within the first few bins of the exponential growth region. The lowest 
energy minima correspond to point defects. The onset of the exponential growth region is 
due do a growing number of defects (mainly of the FFCD type) which lead continuously 
to amorphous structures. The longer simulation in Fig. [10] also shows a second peak at 
around 55 eV. This peak is due to amorphous configurations which are related to a sheared 
crystalline structure. Since we do not relax the simulation cell sheared structures cannot 
relax. 

The reason why we show the MH-DOS over an energy interval which is much larger 
than the interval within which we can obtain a reliable C-DOS is the following. If there 
were good agreement between the C-DOS obtained from different force fields the MH-DOS 
would also agree. As seen from Fig. [11] the MH-DOS obtained from different force fields 
are drastically different and one can therefore conclude that the C-DOS are also drastically 
different. The MH-DOS of all potentials are included in Fig. [TT] Stillinger- Weber and the 
Lenosky Tight-Binding show similar features in the low energy region, e.g., the energy gap 
between the single FFCD and higher excited states and the spike between 7 and 8 eV. While 
the EDIP and Tersoff potential show only a single major peak around 10 eV, both Lenosky 
and Stillinger- Weber have a second peak located at about 35 eV which corresponds, as 
discussed above, to sheared structures. This is due to the fact that for these two potentials 
the C-DOS of unsheared amorphous structures is much lower than for the other potentials 
and the MHM starts therefore sampling higher energy regions corresponding to sheared 
structures. The differences in the C-DOS are responsible for the different speeds with which 
the global minimum is found (see Table llHl) . 

V. CONCLUSIONS 

We have shown that DFT and in particular LDA barrier heights are rather accurate for 
rearrangement processes occurring in silicon clusters. This is good news since the estimation 
of diffusion coefficients and other dynamical properties in silicon systems are frequently based 
on DFT calculations. Since it is well established that DFT schemes give highly accurate 
results for structural properties, i.e., local minima of the potential energy surface, DFT 



19 



calculations are able to provide very reliable potential energy surfaces for silicon systems. 
The bad news is that force fields, that are widely used for dynamical simulations in large 
silicon systems, do not faithfully describe the potential energy surface. With the exception of 
the MEAM based Lenosky force field, all force fields give rise to potential energy surfaces that 
are too rugged. In an extended crystalline environment most force fields greatly overestimate 
the configurational density of states because they give rise to many fake defect structures 
which do not exist in more accurate schemes. The situation is even worse for non-periodic 
systems where the large majority of stable structures are fake. Simulations based on the 
use of a single force field should therefore be viewed with caution and should be verified by 
density functional calculations whenever this is feasible. 
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FIG. 8: (color online) The box plots are based on ten low-lying structures of Sii6- The boxes 
contain the values ranging from the lower to the upper quartiles and the median is represented 
with horizontal lines. The maximum and minimum as well as the mean values are plotted separately 
with crosses and squares respectively. The absolute energy value found after relaxation in DFT are 
represented in (a) with arbitrary origin. Plot (b) shows the difference in DFT energy before and 
after relaxation. Plot (c) shows the average displacements per atom before and after relaxation. 
Unitary transformations of the initial and final structure were performed to diagonalize the moment 
of inertia tensor with respect to each atom. The transformations which resulted in the lowest 
displacement were chosen for the plot. In plots (b) and (c) small values indicate better agreement 
of the potential with DFT results. 
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FIG. 9: (color online) Energies of all potentials along a relaxation path in a DFT calculation 
plotted against the added up norm of the average atomic displacement. The starting configuration 
is a Lenosky force field local minimum. The energies are shifted such that the minimum values 
are set to 0. Lenosky reparame t r ized is an unpublished reparametrized Lenosky MEAM in which a 
single FFCD is stable. 
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FIG. 10: (color online) The minima-hopping density-of-states, MH-DOS, represented by a his- 
togram with 100 bins. 1 x 10 6 and 2 x 10 5 structures were found with the MHM, colored red and 
green respectively. 
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FIG. 11: (color online) The normalized MH-DOS as represented by a histogram consisting of 100 
bins on a logarithmic scale. While 200000 values were used for the classical potentials only 25000 
could be calculated with the Lenosky Tight-Binding scheme. The energy is shifted such that the 
ground state has energy 0. 
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